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ABSTRACT 

We consider the contribution of 3rd and 4th order terms to the power spectrum of 21 cm brightness 
temperature fluctuations during the epoch of reionization (EoR), which arise because the 21 cm 
brightness temperature involves a product of the hydrogenic neutral fraction and the gas density. The 
3rd order terms vanish for Gaussian random fields, and have been previously neglected or ignored. 
We measure these higher order terms from radiative transfer simulations and estimate them using 
cosmological perturbation theory. In our simulated models, the higher order terms are significant: 
neglecting them leads to a > 100% error in 21 cm power spectrum predictions on scales of k > 
1 /iMpc -1 when the neutral fraction is (xh) ~ 0.5. At later stages of reionization, when the ionized 
regions are bigger, the higher order terms impact 21 cm power spectrum predictions on still larger 
scales, while they are less important earlier during reionization. The higher order terms have a simple 
physical interpretation. On small scales they are produced by gravitational mode coupling. Small 
scale structure grows more readily in large-scale overdense regions, but the same regions tend to 
be ionized and hence do not contribute to the 21 cm signal. This acts to suppress the influence 
of non-linear density fluctuations and the small-scale amplitude of the 21 cm power spectrum. In 
alternate models, where the voids are reionized before over-dense regions ('outside-in' reionization), 
the effect should have the opposite sign, and lead to an enhancement in the 21 cm power spectrum. 
These results modify earlier intuition that the 21 cm power spectrum simply traces the density power 
spectrum on scales smaller than that of a typical bubble, and imply that small scale measurements 
contain more information about the nature of the ionizing sources than previously believed. On 
large scales, higher order moments are not directly related to gravity. They are non-zero because 
over-dense regions tend to ionize first and are important in magnitude at late times owing to the 
large fluctuations in the neutral fraction. Finally, we show that 2nd order Lagrangian perturbation 
theory approximately reproduces the statistics of the density field from full numerical simulations 
for all redshifts and scales of interest, including the mode-coupling effects mentioned above. It can, 
therefore, be used in conjunction with semi-analytic models to accurately, and rapidly, explore the 
broad regions of parameter space relevant for future 21 cm surveys. 

Subject headings: cosmology: theory - intergalactic medium - large scale structure of universe 



1. INTRODUCTION 

A frontier in observational cosmology is the detection 
of 21 cm emission from neutral hydrogen gas in the 
high redshift intergalactic medium (IGM) (e.g. Scott 
& Rees 1990, Madau et al. 1997, Zaldarriaga et 
al. 2004; for a review see Furlanetto et al. 2006a). 
These observations promise three-dimensional informa- 
tion regarding the Epoch of Reionization (EoR), con- 
straining the nature of the first luminous objects and 
early structure formation. Indeed, several low-frequency 
radio telescopes (PAST, Pen et al. 20 04; MWA, 
http:/ /web. haystack. mit.edu/arrays/MWA/ Bowman 
et, al. 2006; LOFAR, | htip://www.loiar.org/| and SKA, 
http://www.skatelescope.org/), underway or in the plan- 
ning stages, aim at detecting this signal. 

Detailed theoretical modeling is required to forecast 
constraints, and eventually interpret, the results of these 
observations and to understand their implications for 
early structure formation. The first generation of exper- 
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iments will lack the sensitivity required to make detailed 
21 cm maps, and a statistical detection will be necessary 
(Zaldarriaga et al. 2004, Furlanetto et al. 2004a, Morales 
et al. 2005, McQuinn et al. 2006). One statistic of choice 
is the power spectrum of 21 cm fluctuations, although 
other statistical measures should help in characterizing 
this non-Gaussian signal (e.g. Furlanetto et al. 2004b). 
It will be subtle to infer quantities like the volume-filling 
factor and size distribution of HII regions from the ob- 
served 21 cm power spectrum, and more generally to 
extract information regarding the ionizing sources. 

In this paper, we consider an effect that while impor- 
tant for accurate calculation and interpretation of the 21 
cm power spectrum, has been neglected in many previous 
calculations. The 21 cm brightness temperature involves 
a product of the gas density and the hydrogenic neutral 
fraction, and so the 21 cm power spectrum involves 3rd 
and 4th order terms. In Fourier space these contributions 
involve particular integrals of the 3 and 4-pt functions, 
the bispectra and trispectra of the various fields. We 
will thus loosely call these terms 3 and 4-pt terms even 
though in real space they involve only two points. 

The usual intuition is that on scales much smaller than 
that of a typical HII region, the 21 cm power spectrum 
should be proportional to the density power spectrum. 
This intuition ignores the coupling between large and 
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small scale density fluctuations which arise naturally dur- 
ing the non-linear growth of structure via gravitational 
instability (e.g. Bernardeau et al. 2002). Indeed, mode 
coupling should give rise to non-vanishing 3-pt terms, as 
we detail in later sections of this paper. Qualitatively, 
structure grows more rapidly in regions which are over- 
dense on large scales: an over-dense region acts like a 
closed universe with a boosted matter density. The same 
regions, however, contain more sources and are ionized 
before underdense regions in our models (Sokasian et 
al. 2003, Furlanetto et al. 2004a,c , Iliev et al. 2006, 
Zahn et al. 2006). As the universe reionizes, the large 
scale over-dense regions quickly become 'dark' in a 21 
cm map. One can think of the 21 cm field as a 'masked' 
density field, with the ionization field playing the role of 
the 'mask'. Unlike in a typical galaxy survey, however, 
the mask is itself correlated with the density field, pref- 
erentially removing large scale over-dense regions which 
contain boosted levels of small scale structure. The up- 
shot of this is that, in models where over-dense regions 
are ionized first, mode-couplings should suppress the con- 
tribution of density fluctuations to the 21 cm power spec- 
trum. The aim of the present paper is to demonstrate 
this effect and quantify its importance. 

In 33 we calculate the 21 cm power spectrum (for il- 
lustrative purposes, in real as opposed to redshift space) 
from radiative transfer simulations, and demonstrate the 
significance of 3 and 4-pt terms. In H3.ll we look at the 
small scale effects and motivate the importance of the 
higher order terms with analytic calculations based on 
2nd order cosmological perturbation theory. In H3.2I we 
study the large scale limit of the higher order terms. We 
then illustrate (JU the dependence of the higher order 
terms on the properties of the ionizing sources. Next 
we ($5J refine the fast numerical scheme of Zahn et al. 
(2005) to include the higher-order effects studied here. 
Additionally fffijl- we examine how the results depend on 
redshift and ionization fraction. Here we provide results 
in redshift space, generalizing the illustrative calculations 
of the previous sections. Finally, we discuss our findings 
and conclude in U3 

2. THE 21 CM POWER SPECTRUM 

In this section, we define terms and measure separately 
the contribution of low and high order terms to the 21 
cm power spectrum with radiative transfer simulations. 
For simplicity, we presently neglect the effect of pecu- 
liar velocities which we incorporate subsequently in 
Ignoring peculiar velocities, the 21 cm brightness tem- 
perature, relative to the CMB, at observed frequency, u, 
and redshift, z, is (e.g. Zaldarriaga et al. 2004): 
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In this equation, xh is the hydrogenic neutral fraction, 
1 + 5 P is the gas density in units of the cosmic mean, 
Ts is the spin temperature, and Tcmb is the CMB tem- 
perature. The other symbols have their usual meanings. 
Throughout this work, we will make the usual simplify- 
ing assumption that Ts » Tcmb globally during reion- 
ization, implying ST oc (1 + 5 p )xh, (Ciardi & Madau 



2003, Chen & Miralda-Escude 2003, Furlanetto 2006b, 
Pritchard & Furlanetto 2006). 

In the limit that Ts » Tcmb, and ignoring peculiar 
velocities, the 21 cm power spectrum can be decomposed 
into the sum of several terms (generalizing the formula 
in Furlanetto et al. 2006c): 

Alr(fe) = (T b ) 2 (x H ) 2 [Aljjk) + 2Ai. Sp (k) 

+ 2A! t4pA (*) + 2A^ itfp (A0 

+ Al Wp (k)} (2) 

In this equation 5 X — (xh — {xh})/ (xh) is the fractional 
fluctuation in the hydrogenic neutral fraction, and (Tb) 
is the average 21 cm brightness temperature relative to 
the CMB. Here and throughout A 2 h (k) indicates the di- 
mensionless cross-power spectrum between two random 
fields, a and b. The quantity A^ a (fc) = k 3 P a , a (fc) / '(2tt 2 ) 
is the contribution to the variance of field a per ln(k), 
and this relation defines a dimensional power spectrum, 
-Pa,a(fc), with our Fourier convention. The terms on the 
first line of Equation J2J) are the usual low-order terms, 
representing the power spectrum of neutral hydrogen 
fluctuations, the cross power spectrum between neutral 
hydrogen and gas over-density, and the density power 
spectrum, respectively (e.g. Furlanetto et al. 2006c). 

The terms on the following lines, which we refer to as 
'higher order', are the focus of our present paper. Note 
that the term A 2 p s (k) indicates here the fully non-linear 
density power spectrum. In spite of this, we will loosely 
refer to it, as well as the other terms on the first line 
of Equation as 'low order' since they contain only 
two fields, and because previous analytic calculations in- 
cluded non-linear effects for these terms using the halo 
model (e.g. Furlanetto et al. 2004a). Also note that 
the 4-pt term, A 2 s s s (k), is non- vanishing even in the 
case that S p and 6 X are each Gaussian. In this sense it 
is perhaps misleading to refer to it as 'higher order', but 
we will do so throughout for convenience. 

If we ignore the 3-pt terms, and consider scales much 
smaller than that of a typical bubble, the dominant terms 
from Equation are A 2 - s (fc) and A 2 s s s (k). In the 
Gaussian approximation, we will show in the next sec- 
tion (Equationinj that the term A 2 s s s (k) approaches 

A 2 - s (k) J dlnk 3 A^ x 5 (ka) on small scales; i.e. the ex- 
pression reduces to the variance of the field 5 X multi- 
plied by the density power spectrum. Adding this 4- 
pt term to the density power spectrum piece, we ex- 
pect the small scale limit of Equation @ to approach 
A|i(fc) ~ (T b ) 2 (x 2 H )A 2 s s (k). Further, to the extent 
that Xh is either '1' or '0' at each point in the IGM, 
(x 2 H ) ={x H ), implying A^(fc) oc (x H )A 2 pSp (k) on small 
scales. Note that if we had neglected the 4-pt term en- 
tirely, we would have obtained a different small scale 
limit, A^oc (x H ) 2 A 2 p 5p (k). 

These expressions illustrate the usual intuition that, 
on sufficiently small scales, the 21 cm power spectrum 
should trace the density power spectrum. In fact, we 
will show that the 3-pt terms in Equation @ are gener- 
ally substantial, resulting in large deviations from each 
of these small scale limits. 

2.1. Radiative Transfer Calculations 
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Fig. 1. — Importance of higher-order terms for 21 cm power 
spectrum calculations. Top panel: The black solid line shows the 
simulated 21 cm power spectrum in real (as opposed to redshift) 
space at z = 6.89, at which point the volume- weighted ionization 
fraction is x- l}V = 0.48. The green dotted line shows the contri- 
bution to the 21 cm power spectrum from the low-order terms. 
The magenta dot-dashed line shows the sum of the higher order 
terms for wave-numbers in which the sum is positive, while the 
blue dashed line shows the absolute value of the sum where it is 
negative. Bottom panel: The fractional contribution of the higher 
order terms as a function of scale. On small scales, our results 
differ from the low-order expectation at the ~ 100 — 250% level. 



In order to check the effect of the 'higher order' terms 
we measure each term in Equation (J2J separately, using 
the reionization simulations of Zahn et al. (2006), Mc- 
Quinn et al. (in prep). These simulations follow the 
growth of HII regions in a cubic box of co-moving side 
length L hox = 65.6 /i" 1 Mpc . The Zahn et al. (2006) 
calculations start from an N-body simulation run with 
an enhanced version of Gadget-2 (Springel 2005), track- 
ing 1024 3 dark matter particles, and resolving dark mat- 
ter halos with mass M > 2 x 10 9 M©. Ionizing sources 
are placed in simulated dark matter halos, using a sim- 
ple prescription to connect ionizing luminosity and halo 
mass (Zahn et al. 2006). The simulated dark matter den- 
sity field is then interpolated onto a 512 3 Cartesian grid 4 , 
and we subsequently assume that the gas density closely 
tracks the simulated dark matter density field (see Zahn 
et al. 2006 for a discussion). Finally, radiative transfer 
is treated in a post-processing stage using the code of 
McQuinn et al. (in prep.), a refinement of the Sokasian 
et al. (2001, 2003) code, which in turn uses the adaptive 
ray-tracing scheme of Abel & Wandelt (2002). 

The result of these calculations is shown in Figure 
The black solid line shows the simulated 21 cm 

4 This is higher resolution than the simulations shown in Zahn 
et al. (2006). The higher resolution results will be presented in 
McQuinn et al. (in prep). 



power spectrum in dimensionless units (i.e., it shows 
A2 1 (/c)/(T fc ) 2 from Equation^, the green dotted line 
shows the contribution from the low order terms - i.e., 
it shows the sum of the terms on the first line of Equa- 
tion J5J), and the blue dashed/magenta dot-dashed lines 
indicate the higher order terms. 

The bottom panel further illustrates the importance 
of the 3 and 4 pt terms for accurate predictions of the 
21 cm power spectrum. The curve shows the fractional 
error one makes in neglecting the higher order terms: 
this error is at the ~ 100 — 250% level on scales of 
k ~ 1 — 10 /iMpc -1 . On still smaller scales the simulation 
results are unreliable owing to our limited numerical res- 
olution. As mentioned previously, we might instead have 
included the Gaussian part of the 4-pt function as a 'low 
order' term, in which case the low order terms amount 
to A|i(A;) oc (xh)^s g (k) on small scales. Our results 
differ from this limit by the even larger factor of ~ 250% 
at k ~ 1/iMpc -1 . Clearly, the 3-pt function terms are 
quite important in our simulations. 

3. ANALYTIC ESTIMATES 
3.1. Small scales: perturbative estimates 

In order to develop intuition regarding the 3 and 4-pt 
terms, we will estimate their form analytically, using 2nd 
order Eulerian perturbation theory (e.g. Scoccimarro 
2000, Bernardeau et al. 2002). In this section, we re- 
strict our analytic calculations to small scales (k > lh 
Mpc -1 ), where spatial fluctuations in the density field 
dominate over fluctuations in the neutral fraction, al- 
though we will demonstrate in the next section that the 
higher order terms are generally non- vanishing on larger 
scales as well. In the small scale limit, we can ignore 
mode-couplings between fluctuations in the hydrogenic 
fraction, 6 X , at two different wave-numbers fci and k%, 
which should damp out on scales smaller than the typi- 
cal bubble size. 

First, let us calculate the term Ag^g g (k). Writing 
5 x 5p(k) — ip(k), and using the fact that the Fourier 
transform of a product is a convolution, one has: 

Ps^6 p (k) = {27r) 3 S D (k 1 +k 2 ) (^(k 1 )8 p (k 2 )) 

= 6 D (kt + k 2 ) J ^ (5 B (fc a )« p (fci - fc a )* P (*2)> ( 3 ) 

To lowest non-vanishing order, this expression receives 
contributions from expanding each of the density field 
terms to 2nd order while leaving the remaining ioniza- 
tion field term at 1st order (see Figure for a check on 
the validity of this approximation). First let us expand 
S p (fei — ks) to 2nd order in the linear density field. The 

(2) 

2nd order density field in Fourier space, S p (fci — ks), 
is given by perturbation theory as (e.g. Scoccimarro 
2000): 

-*,)=/ !!||m*i - *» - * - 

xFz( qi ,q 3 )8W(qi)8jP(q a ) (4) 

In this equation, i^Qi, Qa) * s the 2nd order kernel ex- 
pressing mode-coupling from non-linear evolution, Sp 1 ^ 

(2) 

denotes the 1st order density field, and S p denotes the 
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2nd order density field. The second order mode-coupling 
kernel is given by (e.g. Scoccimarro 2000): 

F 2 ( 9l) , 2 ) = | + f^^ + ^ + |(^) 2 (5) 
7 2qiq 2 \q 2 qi J 7 \ qiq 2 J 

Inserting the 2nd order expansion into the integral of 
Equation @, the resulting integrand contains an expec- 
tation value of the form ^S x (k 3 )S p 1 " > (q 1 )S p 1 \q 2 )Sp 1 " > (k 2 ] 

This can be 
uct of two 
each of two 

tributions, 



rewritten as the sum of the prod- 
separate expectation values, with 
terms yielding non-vanishing con- 
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The contribu- 



tion to the integral of Equation ((3J) from expand- 
ing S p (ki — fea) to 2nd order then simplifies to 

25 D (kt + k 2 )P 5p , Sp (k 2 ) J ^F 2 (-k 3 ,-k 2 )Ps x ,s p (k 3 ). 

A second, similar term, arises from expanding S p (k 2 ) to 
2nd order. In total, our expression for this 3-pt term to 
2nd order in perturbation theory becomes 

= 25 D (fei + k 2 )P Sp ,5 p {k 2 ) 

d 3 k 3 
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F 2 {-k 3 ,-k 2 )Ps^ Sp (k 3 ) 
F 2 (-k 3 ,k 3 -k 1 )P Sxt sJk 3 ) 



cP^dfex-feal)- 
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sM dlnk 3 Ai i5 (k 3 ). (7) 



A useful further approximation results from tak- 
ing the small scale limit of this equation. The 
first term simplifies, after spherically-averaging to 
(34/21)P Sp ,s p (ki) f (OnkaA^g (k 3 ). What about the 

2nd term? Since A'j s is small on small scales, the 
small scale limit follows from taking k\ ^> k 3 . In this 
limit one can pull the density power spectrum out of the 
integral. It is important to then perform the average over 
angles, before taking the \k 3 \ — > limit. This eventually 
yields (13/21)P Spt5 Jkx) J dlnk 3 A^ 5p (k 3 ). Summing the 
two terms together, we arrive at a compact expression, 
accurate in the small-scale limit: 
45 

^S x 5 p ,5 p (fcl) ~ -^Ps p 

This approximate expression has an illuminating form. 
The integral over A| x s ^ is a measure of how well the 
ionized regions trace the over-densities. This integral 
can be written as (S X S P ) = — 1 + (xup) / (xh) (p)- Al- 
ternatively, it can be expressed in terms of the volume- 
weighted and mass-weighted ionization fractions, which 
we denote by x; iV and Xi <m respectively, as (S x S p ) — 
(xi iV — Xi. m )/(1 — x- hV ). This term is negative during 
reionization in our simulations: the ionizing sources live 
in over-dense regions and reionize their surroundings, be- 
fore eating out into under-dense regions. Consequently, 
the mass- weighted ionization fraction always exceeds the 
volume- weighted ionization. 

We illustrate this explicitly in Figure [21 (see also Zahn 
et al. 2006), where we show power and cross-power 
spectra for density and neutral hydrogen fluctuations, 
as well as the cross-correlation coefficient between neu- 
tral hydrogen and over-density, which is defined by 
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Fig. 2. — Power spectra of density and neutral hydrogen fluctu- 
ations, and their cross-correlation, in our simulated models. Top 
panel: The red solid line shows the simulated density power spec- 
trum, while the blue dashed line shows the power spectrum of 
fluctuations in the neutral hydrogen distribution. The green dot- 
ted line indicates the absolute value of the cross power spectrum 
between neutral hydrogen fluctuations and density fluctuations on 
scales where it is negative, while the magenta dot-dashed line is 
the same for wave-numbers where it is positive. Bottom panel: 
The cross-correlation coefficient between density fluctuations and 
neutral hydrogen fluctuations as a function of wave-number. All 
results are at x- liV = 0.48 and z = 6.89. The fluctuations in the 
neutral hydrogen distribution are much larger in amplitude than 
the density fluctuations on large scales. The neutral hydrogen fluc- 
tuations are perfectly anti-correlated with density fluctuations on 
large scales, while the two fields become weakly correlated on small 
scales (see text). 



r Sx , Sp (k) = Ps^5 p {k)/[Ps^MPs pl & p (k)] 1 / 2 . Here we 
show results from a simulation output at z = 6.89, 
at which point the simulated volume-weighted ioniza- 
tion fraction is x- hV = 0.48. The figure illustrates that 
on large scales the cross-correlation coefficient is always 
close to —1, a reflection of how tightly the ionized regions 
track over-densities (and hence neutral hydrogen and 
over-density are strongly anti-correlated). On very small 
scales, the cross-correlation coefficient turns slightly pos- 
itive (r ~ 5%), as ionization fronts extend further along 
underdense 'fingers' through the IGM (see McQuinn et 
al. in prep.). Overall, however, reionization is strongly 
'inside-out' in our simulations, as illustrated in the fig- 
ure. Further, note that A^ s ^ and A^ s ^ are much larger 
than the density power spectrum on large scales. This 
implies, from Equation J7J, that the 3-pt function consid- 
ered here will be much larger than the analogous quantity 



for the density field: our effect is generated by A^ s (k) 
on large scales, while the density 3-pt function is driven 
by the large scale density power spectrum, A^ s (/c), 
which is much smaller. 
Indeed, at z = 6.89, the simulated mass-weighted ion- 
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Fig. 3. — 3-point and 4-point terms from our simulations, and 
perturbative calculations. Top panel: The red solid line shows 
the simulated 4-pt function term. The blue dashed line shows the 
2nd order perturbation theory calculation, according to Equation 
l9l . with the simulated density power spectrum as input. The 
difference between the perturbative and simulated results on large 
scales owes to coupling between 5 X modes, an effect ignored in our 
perturbative calculations. Bottom panel: A similar comparison for 
the most important of the 3-pt terms, following Equation J7J. We 
show the absolute value of this term, which is negative over all 
wave-numbers in our simulation. 



ized fraction is 



= 0.61 



giving (8 X 5 P ) the substantial 



value of (S x 5 p ) = —0.25. Furthermore, the pre-factor in 
Equation J7J) is 45/21 and this term enters into the 21 cm 
power spectrum (Equation [21 with an additional factor 
of 2. This demonstrates that 2Ps x s ,s p (k) ~ — Ps„,S p {k) 
on small scales : we expect the 3-pt term to be compa- 
rable in size, and opposite in sign to the density power 
spectrum contribution. 

Strictly speaking, our perturbative expressions contain 
the linear density power spectrum, rather than the fully 
non-linear density power spectrum. However, we expect 
an expression similar to Equation J7J) in the fully non- 
linear regime, containing instead the fully non-linear den- 
sity power spectrum, and with a different, perhaps scale- 
dependent coefficient. Indeed, integrating Equation JJ]J 
numerically with the non-linear density power spectrum 
as input, provides a fairly good estimate of the small 
scale 3-pt term in our simulation. We show this compar- 
ison in the bottom panel of Figure |3] 5 In our fiducial 
model at this redshift, 5 X mode-couplings result only in 
a small correction to the large scale 21 cm power spec- 
trum |QJ , although we will illustrate in 21 an d 3H1 that 
their effect can be significant in other models and at dif- 

5 Note that the pre-factor implied by our numerical integration 
is larger than the 45/21 of Equation J7J, which is accurate only in 
the limit of very small scales. In other words, the 3-pt correction is 
even more important than implied by this approximate expression. 



ferent redshifts. On small scales, the perturbative cal- 
culations provide a fairly good match to the simulation 
results, except they appear to mildly underestimate the 
strength of the 3-pt term at high k. At any rate, our per- 
turbative calculation clearly illustrates that the mode- 
coupling effect will be significant, and demonstrates that 
its strength depends on how well the ionized regions trace 
over-densities. 

Expressions for the other 3-pt term and the 4-pt term 
can be derived in a similar manner. The expression for 
the other 3-pt term is: 
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,(k 3 ). (8) 



This is clearly less important than the above 3-pt term, 
since Ps p ,s x (ki) « Ps p ,s p {ki) on small scales, although 
we find that this term is not completely negligible (it 
contributes at the ~ 10% level on small scales). Finally, 
to lowest non-vanishing order the 4-pt term is given by: 
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The first term will dominate on small scales, and be 
roughly Ps x s p ,s x 8 p (h) ~ Ps p ,s p {k{) J dlnk 3 A^ i<5<B (k 3 ). In 
the top panel of Figure we compare the simulated 4-pt 
term and the results from Equation 0, with the sim- 
ulated Ps s p {k), Ps x ,s x {k), and Ps x ,s p (k) as input. The 
agreement is comparable, although slightly worse, than 
for the 3-pt function. At this redshift, the 4-pt term 
works out to be comparable in magnitude to the domi- 
nant 3-pt term, except that the 3-pt terms come in with 
an additional factor of 2 (Equation [21 and hence repre- 
sent the dominant correction. Further, we have argued 
that this 3-pt correction will be comparable to the contri- 
bution from the density power spectrum, in good agree- 
ment with the simulation results of Figure ^ 

3.2. Large scales 

What about large scales? For the specific case shown 
in FigureE](z = 6.89, x\ jV = 0.48), the higher order terms 
have only a small fractional effect. However, as we detail 
presently and in the higher order terms are gener- 
ally non-negligible even on large scales. In general, it 
is challenging to calculate these terms analytically on 
scales where spatial fluctuations in the neutral hydrogen 
distribution are comparable in strength to density fluc- 
tuations. However, in this section we demonstrate that 
we can analytically estimate the higher order terms at 
high ionization fractions and on large scales, when fluc- 
tuations in the neutral hydrogen distribution strongly 
dominate over density fluctuations. 

In order to gain insight into the properties of the higher 
order terms on large scales, we begin by calculating the 
cross correlation coefficient between 5 x 8 p and each of S p 
and S x . The results of this calculation are shown in Fig- 
ure ^ where we examine a range of redshifts and cor- 
responding ionization fractions. At z — 8.16, when the 
ionization fraction is x- ltV = 0.11 (left panel), the figure 
illustrates that 5 X S P is strongly correlated with fluctua- 
tions in the hydrogenic neutral fraction, S x , and strongly 
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Fig. 4. — Bottom "panel: Cross correlation coefficient between 
fixdp with each of S x and 5 p at (from left to right) z = 8.f6,6.89, 
and 6.56 respectively. The volume- weighted ionized fractions at 
the different redshifts are x\ v = 0.11,0.48, and 0.70 respectively. 
At early and late times, the magnitude of the correlation coeffi- 
cient is large in each case. Top panel: Bias factors, showing the 
ratio of A? s and A? s ; the ratio of A? , , and A? s : and 

o x ,o x Op, Op o x dp,o x o x ,o x 

the ratio A| s s over A? g . At early times (far left panel) and 

on large scales, all of the various power spectra track the density 
power spectrum with scale-independent bias factors. In the inter- 
mediate redshift case (middle panel) a scale independent bias factor 
appears to be a poor approximation. At late times (right panel), 
Ag s g tracks the S x power spectrum closely with a roughly scale 
independent bias factor. 



anti- correlated with the density field, 5 p . Further, ex- 
amining the upper leftmost panel, one can see that the 
power spectrum of neutral hydrogenic fluctuations and 
the 3-pt terms each amount to a scale independent bias 
factor multiplied by the density power spectrum on suf- 
ficiently large scales. At intermediate ionization fraction 
(a^i,v = 0.48, middle panel) the cross-correlation coeffi- 
cient between S X 5 P and each of 6 X , S p is reduced and the 
3-pt terms are correspondingly less important, as already 
illustrated in Figure ^ 

Finally, at high ionization fraction (x iv = 0.70, right 
panel), the correlation coefficients reverse sign and are 
again large in magnitude. The uppermost right panel 
illustrates that the bias between neutral hydrogenic fluc- 
tuations and density fluctuations, ^/Aj^^ is large 
and has a strong scale dependence. The bias between 
5 x 5p and 8 Xl is however, relatively independent of scale. 
This results because, at this ionization fraction, fluctu- 
ations in the hydrogenic neutral fraction strongly dom- 
inate over density fluctuations with A 2 - s /A 2 s ~ 30 

at k ~ O.lh Mpc -1 . In this case S X 5 P just directly tracks 
S x with a roughly scale-independent bias factor. It is 



than the other 3-pt term, A 2 ^ s , on large scales at this 
ionization fraction. It is also significantly larger than the 
4-pt term, A 2 s s s , on the scales of interest, and hence 
represents the dominant higher-order correction at this 
ionization fraction. 

Can we understand analytically the weak scale depen- 
dence and strength of the bias between As x s p j x and 
A5 xi 5 x at late times? The following calculation will be 
facilitated by considering the neutral hydrogenic field x, 
rather than examining fluctuations in the neutral hydro- 
genic field, S x . The relevant 3-pt term, /A 2 xSp X1 is related 
to our usual terms by the equality: 



Jk) = {xhY 



A 2 



(k) . (10) 



We proceed to calculate A xS x (k). It is simpler to 
understand what is going on by thinking about the 
correlation function rather than the power spectrum. 
Thus we are looking for a formula for the correlation 
(x(l)S p (l)x(2)) where 1 and 2 indicate two different 
points. We also recall that we are working in the regime 
where the x fluctuations are much larger than the density 
ones, so correlations will be determined by the structure 
of the x field. For simplicity we will consider the case 
when x can take only the values or 1. We can then 
write, 



(x(l)S p (l)x(2)) = J dS p (l)S p (l) x 

P(S p (l),x(l) = l,x(2) = l) 



P(x(l) = l,x(2) = 1) x 
dS p (l)S p (l)P(S p (l)\x(l) = l,x(2) = 1). 



We now note that P(x(l)=l,x(2)=l) is nothing other 
than (x(l)x(2)). To approximate the integral we make 
use of our assumption that it is the x field that dom- 
inates the correlations so that we can approximate 
P(S p (l)\x(l) = l,a;(2) = 1) « P(S p (l)\x(l) = 1). In 
other words we are neglecting all correlations between 
the density and the neutral fraction at widely separated 
points other than the ones that originate in the correla- 
tions of x. We thus obtain, 



(x(l)6 p (l)x(2)) « (x(l)x(2)} 



(x) ■ 



where we have used the identity: 
dS p S p P(d p \x = 1) 



(x) • 



(11) 



(12) 



Thus, the correlation function (x(l)5 p (l)x(2)) has the 
same shape as that of x, (x(l)x(2)). 

The analytic prediction for the 3-pt term, in the limit 
of strong neutral hydrogenic fluctuations, is hence: 



A„ 



,(k) * (S x 6 p )Al<k). 



(13) 



A very similar argument applies to the 4-pt function, 
yielding: 



A 



(k) ~ (5 x 5 p ) 2 Al x {k). 



(14) 



also clear from the figure that A 2 S s s is much greater 



In Figure we compare the analytic predictions of 
Equations l|13l) and (|14f) with simulation measurements. 
Our simple formulae capture the effect relatively well. 
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Fig. 5. — Analytic predictions for large scale bias at z = 6.56, 
v = 0..70. Bottom panel: Similar to the bottom panel of Figure 
PI this panel shows the cross correlation coefficient between x8 p and 
x. Top panel: The solid black line, the red dotted line, and the 
blue dashed line show various bias factors, similar to the ones in 
Figure[I] The horizontal cyan solid and dashed lines show analytic 
predictions for A 2 xS /A^ and A% g s /A| )a . respectively. 



It becomes clear that in this regime the higher order 
terms have a very different origin: they are not related to 
gravity. They simply result from the large fluctuations 
in the neutral fraction and the fact that in our models 
a neutral point is more likely to be underdense, so that 
(8pX) is non-zero. 

We can also obtain formulae valid in the early regime, 
when it is density fluctuations that dominate the correla- 
tions. The generic expectation value we have to calculate 
to compute the power spectrum involves the probability 
P(S p (l),x(l) = l,6 p (2),x(2) = 1). In the early regime 
we can approximate it by: 



P(6 p (l),x(l) = l,S p (2),x(2) = 1) 
P(x(l) = l\5 p (l))P(x(2) 



P(S p (l),S p (2))x 
1|<U2)), (15) 



which assumes it is the S p correlations that dominate. 
We can then write: 



P(5 p (l), S p (2)) « 



P(S P (1))P(S P (2)) x 



(16) 



with £(r) the density correlation function and a 2 its vari- 
ance. This approximation leads to all power spectra be- 
ing proportional to that of the density. Again this is 
a regime were higher order moments do not depend on 
gravitational non-linearities (as we use only the lowest or- 
der form of the joint distribution function of the density) . 
They are made relevant only by the large fluctuations of 
x. 



Finally, we note that the change in behavior of the cor- 
relation coefficients seen in Figure 0] as reionization pro- 
ceeds can be understood as reflecting the transition be- 
tween fluctuations being dominated by S p at early times 
and being dominated by x at later stages. For exam- 



ple, Equation 1)11) [I shows that A 



S~S D ,S. 



has a contribution 

allu "*«„,x- At early times A^ dominates 
but becomes subdominant late during reionization. This 
transition leads to the change in sign of the cross corre- 

^2 



lation coefficients. Similar arguments apply to A^ 



3.3. Convergence with box-size 

The mode-coupling effects described above imply that 
a rather large simulation volume is needed to accurately 
simulate the small-scale 21 cm power spectrum. In- 
deed, owing to the significant transfer of power from 
large to small scales, one might question whether even 
our small scale 21 cm power spectra results are reliable 
given our limited box-size of L = 65.6 h~ x Mpc . Equa- 
tions 0,®; and © demonstrate that the higher-order 
effects depend on two integrals: J dlnk^A 2 s (k.3), and 

f dinks A 2 s (-^3)1 which represent the cross-correlation 
between S x and S pi and the variance of 5 X , respectively. 

The convergence of our small-scale 21 cm results then 
depend on the convergence of these two integrals with 
increasing box-size. We investigate the convergence of 
(6 2 ) and (5 X 8 P ) using the analytic scheme of Zahn et 
al. (2005, 2006), which provides a quick and accurate 
check. More specifically, we apply the Zahn et al. (2005, 
2006) methodology by generating a Gaussian random re- 
alization of a z = 6.89 density field in a box of side- 
length 250 ft- -1 Mpc and calculate the desired quanti- 
ties. Although using a Gaussian random density field 
ignores the mode-coupling effects of present interest, it 
does yield an accurate calculation of (d 2 ), and (S X 5 P ) 
(Zahn et al. 2006), which determine the convergence 
of our results with increasing box-size. To test conver- 
gence, we compare our calculations of these quantities in 
the 250/i _1 Mpc box to calculations done with smaller- 
volume Gaussian realizations of the same density field. 

The result of this test is shown in Figure for ioniza- 
tion fractions of X{ tV ~ 0.5 at z = 6.89, and x- 1<v ~ 0.7 
at z — 6.56. With our current simulation box-size of 
ibox = 65.6 h~ x Mpc , (S X 5 P ) and (5 2 ) have converged to 
84% and 91%, respectively at Xi yV ~ 0.5, which imply 
that our predictions of the high k 21 cm power spec- 
trum should be relatively free of errors owing to miss- 
ing large scale power. For smaller box-sizes the error 
increases, with a box-size of Lbox ~ 30 h^ 1 Mpc result- 
ing in a ~ 50% error. The sensitivity of (8 X S P ), (S 2 ) to 
large scales arises because HII regions around individ- 
ual, highly-clustered sources quickly merge into 'super- 
bubbles' which become quite large (Sokasian et al. 2003, 
Furlanetto et al. 2004a, c, Zahn et al. 2006). Our point 
here is that this impacts also small scale 21 cm power 
spectrum predictions. Naturally, the convergence prop- 
erties will be worse at higher ionization fractions, when 
the HII regions are typically larger than at our fiducial 
ionization fraction of x^ v ~ 0.5. This is illustrated quan- 
titatively by the bottom set of (green) lines in Figure 

El 

4. DEPENDENCE ON SOURCE PROPERTIES 
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Fig. 6. — Convergence of (8 x Sp) and (S x ) with increasing box- 
size. These terms govern the impact of the higher order terms 
on small scales, and the plot hence gauges the convergence of our 
simulated 21 cm power spectrum measurements on small scales 
with increasing box-size. The blue dotted line indicates the size 
of our current simulation box, Lt, ox = 65.6 h~ 1 Mpc . The red 
lines show the convergence of each term for an ionization fraction 
of Xi v = 0.5, while the green lines show the convergence for x\ v = 
0.7. With our current simulation box-size, we expect that (S X S P ) is 
converged with an accuracy of roughly 90%, while (<5j) is converged 
to ~ 80%. The convergence is a bit worse when the ionization 
fraction reaches x\ v = 0.7. 



The calculations in the previous section provide an- 
alytic understanding regarding the 3-pt functions, and 
demonstrate that their strength depends on how closely 
the ionization field tracks large scale over-densities. In 
this section, we illustrate that this provides additional 
leverage in constraining the nature of the ionizing sources 
and the topology of reionization from 21 cm observations. 
In particular, let us examine the higher order contribu- 
tion to the 21 cm signal for three different models for the 
minimum host halo mass of the ionizing sources. 

Specifically, we consider models in which the ioniz- 
ing luminosity is proportional to the host halo mass, 
with sources residing in halos of mass larger than: i) 
the cooling mass (M min = 1.3 x 1O 8 M at z ~ 7, e.g. 
Barkana & Loeb 2001), ii) A/ min = 2 x 10 9 M Q , and iii) 
M m - ln = 4 x 10 10 A1©, respectively. These are toy mod- 
els, but they span a plausible range of properties given 
our extremely limited observational knowledge regard- 
ing the sources that reionized the IGM. The first model 
assumes that all halos down to the cooling mass (i.e., 
halos with T vlr < 10 4 i"T) contribute to reionization. The 
source prescriptions in ii) and iii) might, for example, 
approximate models in which photo-heating has limited 
the efficiency of star-formation in small mass halos (e.g. 
Thoul & Weinberg 1996, Navarro & Steinmetz 1997, Di- 
jkstra et al. 2004), and diminished their contribution 
to reionization. The formal halo mass resolution of our 
simulation is comparable to the minimum source mass 
in model ii) (Zahn et al. 2006), but in model i) we use 
the results of McQuinn et al. (in prep.), which add lower 
mass sources into the simulation with the appropriate 
statistical properties. 

In Figure we show that the impact of the higher 




k [h Mpc- 1 ] 

Fig. 7. — Dependence of higher order terms on the properties 
of the ionizing sources. Top panel: Similar to the bottom panel 
of Figure we show the fractional contribution of the higher or- 
der terms as a function of scale for each of three different source 
prescriptions. The effect is largest when the IGM is reionized by 
numerous low mass sources. Bottom panel: The cross-correlation 
between over-density and neutral hydrogen fluctuation as a func- 
tion of scale. The strength of the anti-correlation depends strongly 
on the source prescription, and is largest for the cooling mass sim- 
ulation. 



order terms differs significantly for these different models. 
In each case the source efficiency is adjusted to yield 
(x H ) ~ 0.5 at z ~ 7. 6 

From the top panel it is clear that the effect depends 
significantly on the source properties. In the case that 
the minimum mass is M m j n = 4 x 10 10 M Q , the effect is 
only at the ~ 10% level. For sources with a minimum 
mass of M m i n = 2 x 1O 9 M0, as already illustrated in 
Figure HI the effect is at the > 100% level. Finally, with 
sources residing in halos down to the cooling mass, the 
effect is even larger, at the > 200% level. 

The reason for this sensitivity is demonstrated in 
the bottom panel of the figure, which shows the cross- 
correlation coefficient between density and neutral frac- 
tion fluctuations as a function of scale. In models where 
small mass halos contribute to reionization, the ion- 
ized regions track over-densities out to smaller scales. 
There are two reasons for this. First, when the ioniz- 
ing sources reside in more massive halos they are more 
biased and produce larger bubbles at a given ionization 
fraction (Furlanetto et al. 2006c), masking out neighbor- 
ing voids as well as their immediate overdense environs. 
Next, Poisson fluctuations in the abundance of ionizing 

6 More precisely our cooling mass simulation has (xjj) = 0.54, 
our simulation with M m i n = 2 X 10 9 Mq has (xh) = 0.52, and our 



high mass simulation has M n 
all at 2 = 6.89. 



4 x 10 1U M© has (x H ) 



0.53, 
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sources become increasingly important for rare, massive 
halos, limiting the cross-correlation between ionization 
and over-density (Zahn et al. 2006). As a result, the 3- 
pt terms (Equation |7| |5J) become increasingly important 
when the ionizing sources are very abundant and more 
closely trace over-densities. 

We emphasize that each of these models already differs 
in its large scale 21 cm power spectra, owing to differ- 
ences in A 2 ^ 5 x (fc) and Ag ^ K (fc) between these models. 
Our point here is merely that, owing to mode coupling, 
the small scale 21 cm power spectra gives an additional 
handle on distinguishing the different models. This is 
clearly demonstrated by Figure 

These effects might be offset somewhat if recombina- 
tions are more important than in our simulated models. 
Since gas in over-dense regions recombines faster, recom- 
binations act to decrease the tendency for the over-dense 
regions to reionize first. Further, in models where the 
voids reionize first, as might be the case if mini-quasars 
reionize the IGM (Ricotti & Ostriker 2004, although see 
Zhang et al. 2006), the 3-pt terms should be positive 
since the volume- weighted ionization fraction will exceed 
the mass-weighted ionization fraction in these models. 

Finally, notice that the higher order terms contribute 
even on large scales, particularly in the cooling mass 
model, where they amount to a 40% correction (see also 
■ In this model, the strongest large scale contribution 
comes from the Ag SS (k) term which is small and nega- 
tive on small scales, but becomes significant and positive 
near the bubble scale (see H3.2l for comments on the im- 
pact of the higher order terms on large scales). 

5. 2ND ORDER LAGRANGIAN PT AND A NUMERICAL 
REIONIZATION SCHEME. 

In order to forecast the ability of future 21 cm sur- 
veys to constrain reionization physics, we require fast 
and reasonably accurate predictions for the expected 21 
cm signal, spanning a wide range of model parameters. 
For this purpose, rapid semi-analytic calculations are ex- 
tremely valuable. One such semi-analytic scheme is that 
of Zahn et al. (2005), which is essentially a Monte-Carlo 
implementation of the analytic model of Furlanetto et al. 
(2004a). 

The original Zahn et al. (2005) scheme uses Gaussian 
random realizations of cosmological density fields, and 
therefore neglects the mode-coupling effects that are the 
topic of our present paper. This semi-analytic scheme 
can alternatively be applied using the density field, and 
halo distribution, from a full cosmological N-body simu- 
lation, as in Zahn et al. (2006). In this case, the approx- 
imate scheme includes the higher order effects discussed 
above, and the results agree well with more detailed ra- 
diative transfer calculations. Presently, we aim at still 
more rapid calculations, that additionally remove the ex- 
pense of performing an N-body simulation. More specifi- 
cally, we refine the Gaussian random field scheme of Zahn 
et al. (2005) by generating cosmological density fields 
according to 2nd-order Lagrangian perturbation theory 
(2LPT). 

This ambition is plausible given that even relatively 
small scales (k < 10 h Mpc^ 1 ) are still in the quasi- linear 
regime near z ~ 6. 7 Owing to this, we find that parti- 

7 Although non-linear contributions to the density power spec- 



cle distributions set up using 2LPT accurately capture 
the dark matter density distribution at the redshifts and 
scales of interest. 

In 2LPT, as in the Zel'dovich approximation, each par- 
ticle is displaced from its initial Lagrangian position, x, 
to a final Eulerian position, q. The mapping between La- 
grangian and Eulerian positions in 2LPT depends, how- 
ever, on each of the first order potential, 4>^(q), and the 
2nd order potential, (f>^ 2 \q). Specifically, the mapping is 
described by (Scoccimarro 1998): 

x = q - £>i V,^ (q) + D 2 V q ^ (q) (17) 

In this equation, D\ and D 2 are the 1st and 2nd or- 
der growth factors (Scoccimarro 1998). The particle pe- 
culiar velocities satisfy a similar equation (Scoccimarro 
1998). The first order potential <f>^(q) obeys the Pois- 
son equation V 2 ^ 1 ) (q) = S^- 1 ' (q), while the 2nd order po- 
tential satisfies a separate Poisson equation, V 2 ^ 2 ^ (q) — 

E^M^te^teMC] 3 } ( Buchert et aL 1994 > Scoc - 

cimarro 1998). 

The procedure for generating particle realizations of 
a 2LPT density field is then straightforward. First, one 
generates a Gaussian random realization of the first order 
density field with the appropriate linear power spectrum. 
Second, one solves the first order Poisson equation for the 
potential, (f>^ (q). From the first order potential, one 
can calculate the source term in the 2nd order Poisson 
equation and subsequently solve this equation for the 
2nd order potential, 0( 2 )(g). Finally, one displaces each 
particle from its initial cell center using Equation l|17fl 
with the 1st and 2nd order potentials as input. More 
specifically, we generate 2LPT displacements using the 
code of Scoccimarro (1998), Crocce et al. (2006), which 
is now publicly available. 8 

From a 2LPT particle realization, we simply interpo- 
late to find the density field on a Cartesian grid. We can 
then calculate the 21 cm power spectrum, using the ion- 
ization field generated by applying the method of Zahn 
et al. (2005, 2006) to the 2LPT density field, or using 
the simulated ionization field. How well do these predic- 
tions agree with results from our full N-body simulation 
and radiative transfer calculation? We can test this by 
generating a 2LPT realization with precisely the same 
first order displacements as in our N-body simulation. 

Before presenting this comparison, we should make one 
caveat. The method of Zahn et al. (2005, 2006) pro- 
vides a very good, but imperfect match to the ionization 
field from more detailed radiative transfer simulations. 
Presently our aim is to test how well 2LPT predicts the 
higher order terms in the 21 cm power spectrum. We 
thus compare the simulated and 2LPT 21 cm fields, using 
in each case the ionization field calculated from the full 
radiative transfer simulations. We do this comparison, 
rather than using the ionization field from the radiative 
transfer calculation for our simulated 21 cm field, and 
the analytic ionization field for our 2LPT calculation. In 

trum are relatively mild, higher-order contributions to the 21 cm 
power spectrum are substantial, as we demonstrated previously. 
The comparatively larger importance of 3-pt terms for the 21 
cm power spectrum arises because, on large scales, P^ x g (k) 3> 
Ps s p (k). 

http:/ /cosmo. nyu.edu/roman/2LPT/ 
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Fig. 8. — Density power spectra from the 2LPT and N-body 
calculations. Top panel: The red solid line shows the density power 
spectrum from our N-body simulation. The green dashed line is the 
density power spectrum from the 2LPT field. The blue dashed and 
dotted lines show the linear and non-linear density power spectrum 
according to Peacock & Dodds (1996). Bottom panel: The black 
solid line shows the fractional difference between the simulated 
and 2LPT density power spectrum. The red dotted line shows the 
fractional difference between the simulated power spectrum and 
the Peacock & Dodds fitting formula, while the blue dashed line 
shows the fractional difference between the simulated and linear 
theory power spectrum. The green dashed lines indicate fractional 
errors of ±10%, and the cyan dashed line runs through zero. 



Fig. 9. — 21 cm power spectra from 2LPT and N-body calcula- 
tions. Top panel (to be viewed in color): Similar to Figure HI The 
black lines show simulated 21 cm power spectra, the green lines 
shows the contribution to the simulated power spectra from the 
low-order terms, while the blue lines show the absolute value of the 
sum of the higher-order terms. The magenta portions of these lines 
indicate where the sum of the higher order terms is positive, as in 
Figure The solid lines are from the N-body calculation, while 
the dashed lines show 2LPT calculations. Bottom panel: The frac- 
tional difference between the simulated and N-body 21 cm fields. 
The green dashed lines indicate fractional errors of ±10%, and the 
cyan dashed line runs through zero. The agreement is generally 
very good. 



this way, any difference between our 2LPT and N-body 
calculations is attributable to differences in the density 
fields, or the mode-couplings which we study presently. 

First we compare the simulated and 2LPT density 
power spectra. We use 512 3 particles for our 2LPT par- 
ticle realization, drawn from the same initial conditions 
used in our N-body simulation, and interpolate the re- 
sults onto a 512 3 mesh using Cloud-in-Cell (CIC) in- 
terpolation. We have explicitly tested that our results 
are relatively insensitive to the level of shell-crossing in 
this calculation (e.g. Scoccimarro 1998) by filtering our 
initial conditions with several different low-pass filters. 
For our N-body calculation we interpolate our 1024 3 
simulated particles onto a 512 3 mesh. For each of the 
simulated and 2LPT density power spectra, we decon- 
volve the CIC smoothing window before calculating a 
binned, spherically-averaged, density power spectrum, 
and present results on scales only where aliased high- 
k power (e.g. Jing 2005) is insignificant. We do not 
attempt to subtract off shot-noise power from our simu- 
lated results: we have run test simulations with varying 
particle number which indicate that the shot power is sig- 
nificantly sub-Poisson at these redshifts (see also Springcl 
et al. 2005), which makes it difficult to estimate the pre- 
cise level of the shot-noise power. Moreover, even in the 



Poisson case, shot-power would result in a very small 
correction to our 1024 3 particle simulation results : in 
this case, the correction would be only ~ 1% at k ~ lOh 
Mpc -1 , and smaller on larger scales. 

This comparison is shown in Figure 00 demonstrating 
that 2LPT provides a reasonable match to the simulated 
density power spectrum, producing a < 20% underesti- 
mate of the simulated power at k < 5h Mpc -1 and a 
< 40% underestimate at k < 10 h Mpc -1 . As a further 
gauge of the level of agreement, we compare the simu- 
lated density power spectrum with the Peacock & Dodds 
(1996) fitting formula. The Peacock & Dodds (1996) 
power spectrum agrees closely with the 2LPT calcula- 
tion, yet appears to somewhat underestimate the sim- 
ulated density power. We emphasize that the Peacock 
& Dodds (1996) fitting formula is calibrated at z = 0, 
and is by no means guaranteed to match simulations at 
higher redshift. The figure demonstrates that 2LPT is, 
at any rate, just as good as this commonly used fitting 
formula at z ~ 6 and k < lOh Mpc -1 , and is a significant 
improvement over a purely linear calculation. 

Next we compare calculations of the full 21 cm power 
spectrum, separating out contributions from the low- 
order and higher-order terms as in Equation 0. The 
results of this calculation are shown in Figure [5] demon- 
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strating that the 21 cm power spectrum calculations 
agree well (to better than 10% at k < lOh Mpc^ 1 ). The 
agreement is hence even better than the agreement found 
in our density power spectrum calculations (Figure [SJ) . 
This superior agreement, however, owes somewhat to a 
cancellation which arises because the low-order terms in 
our 2LPT calculation are smaller than in our N-body cal- 
culation, while the magnitude of the (negative) higher- 
order terms is also smaller in the 2LPT calculation, as 
illustrated in Figure The agreement might, therefore, 
be a little worse at different ionization fractions, or for 
different models, where this close cancellation may not 
occur. On still smaller scales, where one is dominated 
by the 1-halo term in the density power spectrum (e.g. 
Cooray & Sheth 2002), 2LPT will break down further. 
Additionally, our assumption that gas closely traces dark 
matter should break down (Gnedin & Hui 1998). How- 
ever, even next generation 21 cm experiments such as 
SKA will likely be limited to scales of k < 10/iMpc -1 
(McQuinn et al. 2006, Bowman et al. 2006). In con- 
clusion, 2LPT provides a significant improvement over 
the Gaussian random field scheme of Zahn et al. (2005), 
while requiring very little additional computational over- 
head. 

6. REDSHIFT EVOLUTION AND RESULTS IN REDSHIFT 

SPACE 

In the previous sections, we characterized the impact 
of the higher order terms, ignoring the effect of peculiar 
velocities, and focusing on a single redshift for simplicity. 
In this section, we expand on these calculations by incor- 
porating the effect of peculiar velocities (e.g. McQuinn 
et al. 2006), and by examining the dependence of our 
results on ionization fraction. 

In each case we calculate the spherically averaged 21 
cm power spectrum; i.e. the redshift space monopole, 
as in Zahn et al. (2006). In the case of the monopole, 
the sum of the low-order 21 cm power spectrum terms 
is given by A| lilow = <T b ) 2 (^[A^ + (8/3) + 

(28/15)A^ s } (McQuinn et al. 2006, Zahn et al. 2006). 
In principle, one can write down an analog to Equation 
in redshift space, incorporating all of the relevant 
higher order terms. In practice, we instead simply cal- 
culate the full redshift space 21 cm power spectrum, and 
compare with the decomposition above. 

The results of this calculation are shown in Figure ITU1 
for a range of different redshifts and ionization fractions. 
In the top panel we show the redshift space monopole in 
units of (mK) 2 , which can be more easily compared with 
observational noise estimates than the dimensionless 21 
cm power, which we plotted previously for simplicity. 
The top panel shows the usual qualitative behavior for 
the 21 cm power spectrum redshift evolution (Furlanetto 
et al. 2004a, Zahn et al. 2006): spatial fluctuations in 
the hydrogenic neutral fraction imprint a knee in the 21 
cm power spectrum on large scales when the HII regions 
become sufficiently large. 

The bottom panel of the Figure shows the fractional 
effect of the higher order terms on the redshift space 
monopole. First, let us focus on the small scale behav- 
ior. Focusing on the z = 6.89 curve for the moment, 
and comparing with its real space analogue (Figure ^) , 
it seems that the effect is slightly enhanced in redshift 
space, resulting in up to a factor of ~ 3 difference with 
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Fig. 10. — Redshift evolution of the 21 cm power spectrum 
monopole in redshift space, and fractional importance of higher 
order terms. Top panel: Spherically averaged 21 cm power red- 
shift space power spectrum for different redshifts and ionization 
fractions. Bottom panel: Fractional importance of the higher or- 
der terms for the spherically averaged 21 cm redshift space power 
spectrum. 



the low-order calculation. This presumably results be- 
cause peculiar velocities introduce additional 3-pt terms, 
above and beyond the ones in Equation (J2J , that are ne- 
glected in the usual decomposition. 

At high redshift, z = 8.16, when the volume- weighted 
ionization fraction is x\. v = 0.11 in our model, the small- 
scale suppression is less significant than at z = 6.89. This 
results because the (5 X 8 P ) is smaller at this redshift than 
at our fiducial redshift ((5 X S P ) — —0.12, compared to 
(8 x Sp) — —0.25), amounting to a factor of ~ 2 less sup- 
pression (see Equation[7J) . Finally, at the highest redshift 
considered here, z = 6.56, where the volume- weighted 
ionization fraction is Xi jV — 0.70, the small-scale sup- 
pression is less significant again. This occurs because 
the 4-pt function becomes more significant at low neu- 
tral fraction, and partly compensates the 3-pt function 
suppression. Roughly speaking, the 4-pt function term is 
generated by (S 2 ,) while the 3-pt function term is gener- 
ated by (S X S P ) (see Equations |5] and |7J). Since (5 X ) grows 
more rapidly with decreasing neutral fraction than (5 x S p ) 
(Zahn et al. 2006), the 4-pt function becomes more sig- 
nificant relative to the 3-pt term at x- ljV > 0.5, reducing 
the small-scale suppression, as seen in the figure. 

What about the effect of the higher order terms on 
large scales? The figure clearly shows that the higher 
order terms are generally non-vanishing even on quite 
large scales, with the smallest effect occurring at our fidu- 
cial redshift of z = 6.89, when the ionization fraction is 
%i.v ~ 0.5. Second, note that the dependence on ion- 
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ization fraction/redshift is not monotonic, with low ion- 
ization fractions resulting in a suppression of the 21 cm 
power spectrum signal, intermediate ionization fractions 
yielding an enhanced signal, and the high ionization frac- 
tion leading to a suppressed signal on large scales again. 
These large scale effects are unrelated to gravitational 
instability and have their origin in the fact that in our 
models high density regions tend to ionize first. Just as 
in the small scale case, their amplitude is relatively large 
because the fluctuations in the neutral fraction are large. 
The behavior of the higher order terms on large scales 
depends on whether it is the density or x that dominates 
the correlations. The cross terms have different signs in 
both limits so they go through a minimum around the 
mid point of reionization (see i|HJl . 

7. CONCLUSIONS 

In this paper, we demonstrated the significant impact 
of higher order terms on 21 cm power spectrum predic- 
tions. While on small scales the effect originates in the 
mode-mode coupling induced by gravitational instability, 
on large scales it is unrelated to gravity. It originates in 
the fact that high density regions tend to ionize first. We 
showed that these effects can help distinguish between 
different models for the ionizing sources, and constrain 
how well ionized regions trace large scale over-densities. 
Finally, we demonstrated that these effects can be cap- 
tured in semi-analytic calculations by using 2nd-order 
Lagrangian perturbation theory. 

How important are these effects for upcoming 21 cm 
surveys? The first generation of 21 cm experiments, such 
as MWA and LOFAR, will likely be sensitive only to 
modes with k < l/iMpc -1 (Bowman et al. 2006, Mc- 
Quinn et al. 2006) so only the large scale effects we 
discussed will be relevant. We showed that the higher 
order terms can be significant on larger scales - to the 
extent that ionized regions trace large scale overdensities 



- providing important information regarding the mor- 
phology of reionization. The subsequent generation of 
experiments, like SKA, should extend power spectrum 
measurements to k ~ 10/iMpc -1 (Bowman et al. 2006, 
McQuinn et al. 2006), in which case our small scale ef- 
fect should be extremely important, unless the IGM is 
reionized by very rare sources. 

Another natural question is: is there a statistic that 
isolates the mode-coupling effect described in this paper? 
As the 21 cm power spectrum is produced by several un- 
knowns simultaneously, such a statistic would help iso- 
late the degree of correlation between the density and 
ionization fields. We examined higher-order statistics in 
the vein of Zaldarriaga et al. (2001) to try and isolate our 
effect, but we find that the 21 cm 3-pt function involves 
terms oc {P6^sAk)\ aIge (P Sf>tSp {k)) amall - where 'large' 
and 'small' refer to large and small scales respectively 

- which are non-vanishing even in the absence of corre- 
lations between the density and ionization fields. This 
makes it difficult to disentangle the effect of interest, but 
further work in this direction might be interesting. 

In future work, we will forecast constraints for upcom- 
ing 21 cm surveys using the 2LPT scheme presented here 
(Zahn et al. in prep). Finally, it might be interesting to 
investigate whether analogous 3-pt terms are important 
for accurate predictions of the kinetic Sunyaev-Zel'dovich 
effect. 
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